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We report the results of a Monte Carlo study of a model of (III,Mn)V diluted magnetic semi- 
conductors which uses an impurity band description of carriers coupled to localized Mn spins and 
is applicable for carrier densities below and around the metal-insulator transition. In agreement 
with mean field studies, we find a transition to a ferromagnetic phase at low temperatures. We 
compare our results for the magnetic properties with the mean field approximation, as well as with 
experiments, and find favorable qualitative agreement with the latter. The local Mn magnetization 
below the Curie temperature is found to be spatially inhomogeneous, and strongly correlated with 
the local carrier charge density at the Mn sites. The model contains fermions and classical spins and 
hence we introduce a perturbative Monte Carlo scheme to increase the speed of our simulations. 



PACS: 75.50.Pp, 02.70.Uu 
I. INTRODUCTION 

III-V diluted magnetic semiconductors (DMS) have 
become a very active area of research due ±a their in- 
teresting magnetic and transport properties .lira Thus far 
Gai-^Mn^As has received the greatest amount of atten- 
tion, due to the observation that it becomes ferromag- 
netic withj-a. Curie temperature as high as 110K, when 
x = 0.053.u ! Q More recently, ferromagnetism aJaeve room 
temperature has been observed in (Ga,Mn)N.En3 

The observed ferromagnetism is widely accepted to be 
due to a charge carrier mediated coupling between the 
Mn spins. Several models have been proposed to explain 
the detailed phenomenology of these compounds .LTcJ In 
particular, we have initiated an effort to understand the 
effects of disorder- in—Mai. positions on the properties of 
these compounds. Q'EfEjLJ An additional source of disor- 
der is the large degree of compensation seen in these ma^. 
terials, which has been attributed to As antisite defects,E£l 
because of which, the carrier density is significantly less 
than the Mn density. 

To this end, we study the low Mn concentration regime, 
near and below the metal-insulator transition (MIT) at 
x ~ 0.035, where disorder effects would be expected to be 
the most pronounced. In this limit, we model the charge 
carriers in terms of an impurity band comprised of states 
around each Mn acceptor, which is taken to be the source 
of charge carriers mediating ferromagnetism. Evidence 
for the relevance of an impurity band has been provided 
by a number of recent experiments, such as an scanning 
tunneling microscope studyEa which showed the existence 
of an impurity band in (Ga,Mn)As samples with x = 
0.005 — 0.06. Angle resolved photoemission spectroscopy 
in a sample with x = 0.035 has also revealed a well formed 
impurity band—and confirmed that the Fermi energy lies 
in its vicinity.EB 

A complete, detailed description of the impurity band 
in the presence of compensating defects is an extremely 



difficult enterprise] 3 !] hence we have proposed a simple ef- 
fective Hamiltonian that captures (at least qualitatively) 
the relevant impurity band physics. This Hamiltonian 
was studied comprehensively at the mean-field level, and 
the magnetic propertiesjwere observed to have a number 
of surprising properties. Lru The mean-field magnetization 
curves have very unusual, concave upward shapes, unlike 
the magnetization curves of conventional fcrromagnets. 
Some of these features have been seen in experimental 
measurements, especially for samples with, a low carrier 
density and high degree of compensation!!!! This unusual 
shape of the magnetization curves was identified to be a 
direct effect of positional disorder of the Mn ionsB 

The mean field calculation suggested that there is con- 
siderable inhomogeneity in the magnetization of individ- 
ual Mn spins at temperatures below T c , particularly for 
small values of x. Experimentally, disorder appears to be 
relevant even in the metallic phase - Barkhausen jumps 
were observed in a sample with x = 0.047,E3 indicating 
the presence of frozen-in magnetic disorder. Theoreti- 
cally, disorder has also been credited as leading to an 
instability, in collinear Mn ground states in the metallic 
phase fjEi via spin waves, which can lead to a reduction 
in the saturation magnetization at low temperatures as 
observed in annealed samples with x = 0.05S3 

Increased disorder was also shown to lead to an in- 
creased value of the critical temperature T c in the mean 
field study.ll It is well known that the mean-field approx- 
imation underestimates the effect of thermal fluctuations 
and therefore overestimates the value of the critical tem- 
perature. As a result, it is important to check the mean 
field results against Monte Carlo simulations which prop- 
erly account for thermal fluctuations, to test to what ex- 
tent the phenomenology found in the mean-field study is 
maintained when these fluctuations are included. 

In this paper, we report the results of such a Monte 
Carlo study. We find that the magnetization curves re- 
tain their unusual shape even when thermal fluctuations 
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are included and that, as expected, the critical temper- 
ature is lowered from the mean field value. The extent 
of this decrease is most pronounced for low values of x 
and high compensation. The observation that disorder 
can lead to a. higher critical temperature than a purely 
ordered caseEl is also confirmed, although the change is 
more modest than in the mean field study. The spa=j 
tial inhomogeneity observed in the mean field studiesEl 
is confirmed and found to be little changed by thermal 
fluctuations. A correlation between regions of larger lo- 
cal charge density and larger local magnetization is also 
established. 

This paper is organized in the following manner: in 
section |0] we describe the model Hamiltonian and the 
values chosen for various parameters of the system. We 
also discuss the use of finite size scaling to determine the 
ferromagnetic transition and describe the quantities that 
we calculate. In section [II, we introduce a perturba- 



tive scheme to speed up the Monte Carlo simulations, 
describe our implementation and discuss testing on a toy 
model for which we can compare with exact results. We 
present the results of our Monte Carlo simulations for 
the impurity band model of DMS in sections |[y] and [v|. 
Section ^ summarizes our conclusions and discusses the 
implications of our results for experiments and modelling 
of III-V DMS. 



II. MODEL 

A. Hamiltonian 

When manganese is introduced into GaAs, the Mn im- 
purities have been shown to substitute on the Ga FCC 
sub-lattice of the zinc-blende structurc.pf the undoped 
semiconductor for small values of £.00 However, at 
larger values of x (> 0.07), Mn ions pan-iprm MnAs clus- 
ters, which have the NiAs structure£aE3 Based on these 
experimental findings, we assume a zinc-blende structure 
and Mn substituting only on Ga sites for the low Mn 
concetrations we study. Each Mn impurity has a spin-| 
from its half-filled 3d shell. The nominal valence-II of 
Mn implies that when it substitutes for the valence-III 
Ga, it acts as an acceptor. Thus an isolated Mn can 
bind a hole in an impurity level that we characterize by 
a hydrogenic orbital with Bohr radius as, with wave- 
function 4>i(r) ~ exp[— |r — Ri|/as]. In Mn dgped GaAs 
there are substantial central cell correctionsEZI which we 
phenomenologically incorporate by adjusting the Bohr 
radius a^. 

Whilst the true carriers in this system are holes with 
spin |, we consider the case of electron doping, so the car- 
rier spin is i. This leads to some differencesliij (in par- 
ticulars-trie frustration effects recently claimed for hole 
dopingo are absent in our model). However, the effects 
of disorder and impurity potentials, which are oujj-.main 
focus here, should not be qualitatively changcdflu The 



carriers (electrons) are then in an impurity band below 
the conduction band minimum of the host semiconduc- 
tor. 

We |-jstudy the model introduced by Berciu and 
BhattEQ for which the Hamiltonian is 
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The random positions of the Mn impurities are labelled 
by Ri, and the spin of the Mn impurity at R^ is S^. 
In the electron formalism for the charge carriers, cJ CT 
is the creation operator for a carrier with spin er in 
the bound state associated with the i th impurity. The 
hopping matrix is given by tu = i(|R-i — Rj|)j where 
t(r) = 2(l + r/a B )e~ r / aB Ry0and the Rydberg (Ry) is 
the binding energy of a hole. We assume that the Mn 
spins are strongly localized and hence the exchange inte- 
gral is given by J io = Jo\(f>ij\ = Jo e ~ 2 l Rl ~ R J l/ aB , which 
is proportional to the charge density at the i th Mn site 
of the carrier in the hydrogenic wavefunction around the 
j th Mn site. The external magnetic field is given by H, 
the Lande g-factors for the Mn spins and carriers are g 
and g* respectively, and is the Bohr magneton. Our 
simulations are for zero magnetic field. 

We consider finite size cubic samples with L cubic unit 
cells of the GaAs structure per side, with periodic bound- 
ary conditions, and Nd Mn impurities. Thus x is given by 
Nd/(4:L 3 ), and is related to the concentration of Mn im- 
purities, riMn through riMn — 4x/a 3 , where a = 5.65 A 
is the GaAs lattice constant. The number of carriers 
is Nh — pNd, with p between 0.1 and 0.3, as indicated 
by experimental .studies of samples grown by molecular 
beam epitaxyEJo As mentioned previously, these low 
values of p are due to compensation processes, in ■\sthich 
As antisites are believed to play an important rolc.Ej 

The only difference between our model and that in Rcf. 
f?]|| is that instead of studying spin-| Mn spins, we treat 
the Mn spins as classical Heisenberg vector spins. This 
should be a reasonable approximation since S = 5/2 is a 
large spin and a Quantum Monte Carlo calculation would 
not gain much due to the uncertainties in the materials 
parameters. [At the mean field level, this approximation 
has the effect of lowering T c by a factor of I + -| relative 
to the T c for quantum spins] . One quantity which differs 
significantly wilh-classical spins is the specific heat at low 
temperaturesjfjed which we have not studied. 

The values we use for numerical parameters are as — 
7.8 A, I Ry = 112.4 meV, and Jo = 15 meV, as discussed 
in Ref. ||. With these parameters we have an impurity 
band whose bottom lies around 200-300 meV (2-3 Ryd) 
below the host conduction band and the Fermi energy 
varies from around 13 to 55 meV (depending on x and 
p) above the bottom of the impurity band for the pa- 
rameter range considered. These values are of the same 
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order of magnitude as the observed splitting and band- 
width of the impurity band in angle resolved photoemis- 
sion experiments.E3 

Our goal in this work is partly to understand the de- 
viations from mean- field theory in the model Eq. (^), 
therefore we also perform mean field calculations using 
Langevin functions to represent the polarization of clas- 
sical spins (in the quantum case one uses Brillouin func- 
tions), with which we compare our Monte Carlo results. 



B. Method of Simulation 

Consider a system with classical spin degrees of free- 
dom and fermionic degrees of freedom, such as described 
by the Hamiltonian (|l|) . The assumption of classical spins 
means that one can parameterize the spin at each site by 
its z-component and azimuthal angle, i.e. S, = (Sf, <fii), 
and 



Sf = S\ 1 



st 
s 2 



st 
s 2 



(2) 
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For any given configuration of classical spins, {S} = 
{5 Z ,0}, the Hamiltonian (0) can be diagonalized to give 



H({S, al a n }) = J2 E n ({S z ,cj>})ata n , 

n 

where the states n are a diagonal basis and 



(4) 
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are linear combinations of the c ia and Ci a operators. Here 
ipnaii) is the eigenfunction for the n th energy level at site 
i for spin cr. Using the grand-canonical ensemble, the 
partition function is 



'N d „i 

n/ ds i 



E - E 

ni=0,l njy — 0,1 
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where — 0, 1 is the occupation number of level k (there 
are N = 2Nd levels), and f3 — l/(fcsT) is the inverse 
temperature. Summing over fermion degrees of freedom 
leads to the result 
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This can be cast in a similar form to that used for a 
spin-only system 
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n 



dSf 
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where the carrier free energy T C {{S Z , (/>}) for a given con- 
figuration of classical spins is 



Nu 



F C ({S*, *}) = ~~ E lo g (l + e-^"^.*"-")) • (9) 



The chemical potential fi is determined from the condi- 
tion 



(N c ) = 



d 



where 
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is the Fermi distribution function and (N c ) is the expec- 
tation value for the number of carriers. 

Equations (0) and (^) imply that we can use Monte 
Carlo techniques to evaluate various thermodynamic 
quantities, except that we must use the carrier free en- 
ergy rather than the internal energy in the Metropolis 
algorithm. An important point to note is that since the 
calculation is in the grand canonical ensemble, both the 
temperature and chemical potential must be kept fixed 
during the simulation. This differs from a recent Monte 
Carlo studyu3 in which N c was held constant at each 
Monte Carlo step by varying fj, rather than holding it 
constant, in our simulations we hold f3 and [i fixed for 
each run and then we average over disorder using sam- 
ples with equal (N c ), since we wish to average samples 
which have the same x and p. 



C. Magnetic and Thermodynamic Quantities 

Our Monte Carlo simulations allow us to compute vari- 
ous magnetic and thermodynamic quantities after equili- 
bration. We perform equilibrium thermal averages (indi- 
cated by (...)) for each sample (realization of disordered 
Hamiltonian) and then average over many different real- 
izations of disorder. (The disorder average is indicated 
by an overbar 7TT ) . We collect data for both global (bulk) 
and local quantities. We first consider the global quanti- 
ties we study. We compute the moments of the average 
magnetization per Mn spin, M q , and the average magne- 
tization per carrier, m q , at each Monte Carlo step, where 
the q th moments are given by 
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In the above equation S = § . The average spin per Mn, 
Smui an d the average fermion spin, s c are thus 

S Mn = W), sc = -(m). (14) 

Note that with our definition, Smti is normalized, i.e. for 
a fully polarized spin state Smu — L The negative sign 
for the fermion magnetization is due to the antiferromag- 
netic interactions between Mn spins and carriers which 
leads to oppositely oriented polarizations. The Mn sus- 
ceptibility is given by 
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(m 2 ) - (My 



while the carrier susceptibility is 



Xh = (3 



(m 2 ) — (mY 



(15) 
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For each measured quantity, we calculate the statistical 
errors in the conventional manner from the variance. One 
of the limitations of studying a magnetic model on small 
lattices, as we are forced to do here, is that it is not easy 
to identify the position of a thermodynamic transition 
from considering quantities such as the magnetization. 
However, from finite size scaling theory, quantities that 
are dimensionless are particularly useful to identify the 
transition temperature. One such quantity, known as the 
Binder cumulant, is given by 



G(L,T) = \ 
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(M 4 ) 
(M 2 ) 2 
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G(L, T), which measures the ratio of the fourth moment 
to the square of the second moment of the magnetization 
for a finite system, is defined such that in the param- 
agnetic phase G(L,T) decreases with L, and tends to 
zero as L — » oo, while in the ferromagnetic phase it 
increases with increasing size L and tends to unity in 
the thermodynamic limit. Near the transition tempera- 
ture T c , being dimensionless, it has th t ; fi aite size scal- 
ing form G(L, T) = G[L»(T - T c )],BEl43 where v is 
the exponent of the diverging spin-spin correlation length 
£ - (T-T c y v . Consequently, at T c , G(L, T c ) is indepen- 
dent of L; T c can be identified by a simultaneous crossing 
of G(L, T) vs. T curves for different L. Because G is di- 
mensionless, and depends only on the ratio L/£, rather 
than both L and £, this method is found to be more 
reliable in determining T c than analysis of the onset of 



magnetization, or peaks in the magnetic susceptibility in 
finite sized samples. 

Local quantities we calculate are the local charge den- 
sity at each Mn site: 



Pi 



3,° 



(18) 



and the local magnetization, which we define to be the 
average projection of the spin at site i on the total mag- 
netization M 



M local = (g. . M ^ 



(19) 



We are interested in the individual distributions V(pi) 



local \ 



to 



and V(M\ ocaX ) and the joint distribution Vj(p l , M, 
characterize the local environment at different Mn sites. 



III. PERTURB ATIVE MONTE CARLO (PMC) 

In principle, to perform Monte Carlo simulations on a 
model with fermion and classical degrees of freedom, one 
needs to diagonalize the fermion part of the Hamiltonian 
after each "spin flip" (more precisely a spin rotation), i.e. 
a new choice of the classical variable. This leads to new 
eigenvalues which are used to compute the change in the 
carrier free energy J- C ({S Z , </>}). This is computationally 
time consuming, and hence one would prefer a quicker, 
approximate method which is still reasonably accurate. 
One such approach is the hybrid Monte CapLp (HMC) 
algorithm used recently on the DMS problem,E3 and also 
on the double exchange model.o We have developed an 
algorithm that works in a similar spirit to HMC, that we 
describe below. 

We select a spin, at site i and allow it to perform a 
small rotation: 



St 



SI + 5S? 



(20) 
(21) 



where SSf E [— A/2, A/2], and Sfy S [— \ir, \ir], are re- 
stricted to a small region on the surface of the unit sphere 
(A 1), in a way that leads to uniform sampling. We 
use perturbation theory to compute the change in the 
carrier eigenenergies due to the rotation of the spin at 
site i: 



E n ({S z ,^)^E n ({S z ,^) + 5E n 



(22) 



where 



6E n = -5Si ■ ^ JtJ^^aU^aplpnfiU)- (23) 



a/3 



As a result, SE n ~ O(X) and therefore perturbation the- 
ory is essentially exact as A — > 0. 

We now use Eq. @ to compute the change in carrier 
free energy associated with the spin rotation, and the 
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Metropolis criterion to decide whether to accept the spin 
rotation. We perform such an update for each spin in the 
system using the perturbation scheme. After a complete 
sweep through the system, we compute the eigenenergies 
E n and eigenfunctions ip n (i) corresponding to the new 
spin configuration using exact diagonalization, and start 
a new Monte Carlo sweep. We found that this approach 
was quicker than if we diagonalized after every spin flip, 
generally by a factor of 3-4. 

A. The chemical potential 

One important issue is the choice of the chemical po- 
tential /i for the desired average number of charge carri- 
ers, (N c ) — Nh- To determine the chemical potential we 
consider two replicas of the system, one starting from a 
fully polarized (ferromagnetic) configuration, the other a 
purely random (paramagnetic) configuration. After ev- 
ery few Monte Carlo steps (in practice 5 MC steps worked 
well), we use the condition 

^=£i^Wi' (24) 

n 

to update the value of \x for each replica. When the mag- 
netization M and chemical potential of each replica agree 
to within 2 %, we continue to calculate the chemical po- 
tential for each replica after every 5 Monte Carlo steps, 
but their dynamics are determined by the mean chemi- 
cal potential, /2. This average chemical potential is free 
to vary up until some equilibration time, and then the 
chemical potential is fixed as /i = (ft) where the aver- 
age is over the time during the equilibration period for 
which ft is used to calculate the carrier free energy. The 
equilibration time used depended on the temperature, x 
and p, but was generally between 20000 and 40000 Monte 
Carlo steps. (The magnetization generally equilibrated 
within 2000 Monte Carlo steps, while the remainder of 
the equilibration was required to obtain an accurate value 
of the chemical potential). The fixed chemical potential 
is used for the remainder of the run, during which data is 
collected. We typically use 20000 to 40000 Monte Carlo 
steps to collect data. This procedure was found to ob- 
tain a chemical potential that yielded values of (N c ) for 
the sample that were generally within 2% of Nh (not 
surprisingly it was found to be more effective in larger 
samples where the relative sizes of the fluctuations are 
correspondingly smaller). 

In several Oth er ., st udies of models with fermions cou- 
pled to spinal3't3c3 the electron occupation numbers 
have been taken to be those corresponding to T = 
(i.e. the Fermi distribution is replaced by a Heaviside 
distribution). In our simulations, we allow the filling to 
change as a function of temperature (i.e. we do not as- 
sume degenerate electrons). We do, however truncate the 
number of states we include - their number varies as a 
function of temperature, such that states of high energy 



which have 10~ 5 or lower filling are discarded (the cutoff 
is sufficiently small that the results are not sensitive to 
its value). 

B. Testing the PMC algorithm 

We tested our perturbative Monte Carlo algorithm on 
a simple model of fermions coupled to classical spins 
which admits an exact solution. The Hamiltonian, shown 
below in Eq. ( p5| ) , corresponds to fermions hopping along 
a 1-dimensional chain with TV sites and periodic bound- 
ary conditions. At every lattice site i, there is a classical 
spin Si , and the fermions are equally strongly coupled to 
all spins on the chain. The model Hamiltonian is thus 

n = t J2 c i<? c j° + Jj J2 Sl ' J2 c j<*\ a <*P c JP> (25) 

where energy is measured in units of t = 1 and only 
nearest neighbour hopping is allowed. Note that the ex- 
change scales as J/N to ensure an extensive energy. The 
eigenvalues for this model can be calculated exactly as 
discussed in Appendix Exact results for any num- 
ber of classical spins, N, and Nh fermions can be cal- 
culated, particularly for the magnetization per spin, S, 
the Binder cumulant G(L) and the fermion magnetiza- 
tion s c . A comparison of the exact results for S(T) and 
those calculated using PMC are shown in Fig. [I] for t = J, 
N = 20, N h = 3 and using 40000 Monte Carlo steps. 




FIG. 1. Comparison of exact result (line) and Monte Carlo 
(points) for the magnetization of the toy model, for t = J, 
N — 20 and N h = 3. 40000 Monte Carlo steps were used to 
generate the MC results. 

In Fig. |^ we plot the relative error in the fermion mag- 
netization against the number of Monte Carlo steps (t) 
for the same parameters used for Fig. [I] at T j J = 0.0064 
where the magnetization is about w 0.5. At long times 
the convergence is of the form t~ 2 , where t is the number 
of Monte Carlo steps, as expected for Monte Carlo sim- 
ulations. We have also made a comparison between the 
PMC method and diagonalization after every spin rota- 
tion for the model of DMS that we are really interested 
in [Eq. ([!])], which gave agreement to within the 2 % 
statistical error bars. There is also some contribution to 
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the systematic error from our approach, but we choose A 
to be small (generally A < 0.03), so that this error has 
little effect on our results. 




10000 le-t<)5 
No. of Monte Carlo steps 



FIG. 2. Convergence of the fermion magnetization for the 
toy model with the same parameters as in Fig. |lL at termpera- 
ture T/J = 0.0064. The straight line has a dependence, 
where t is the number of Monte Carlo steps. 



IV. RESULTS: GLOBAL QUANTITIES 

Four different values of x and p were chosen for the 
Monte Carlo study, representative of the concentrations 
and compensation seen in the experimental materials. 
These were: x = 0.01, p = 0.1; x = 0.01, p = 0.3; 
x = 0.03, p = 0.1 and x — 0.03, p — 0.3. Note that using 
Mott's criterion as in Ref. ||, the x = 0.01, p = 0.1 sam- 
ple is in the insulating phase, whilst the remainder of the 
cases considered are somewhere on the metallic side of the 
MIT. In the case of the most metallic sample x = 0.03, 
p = 0.3 it may be inappropriate to neglect the host band 
states for a quantitative description of the magnetic prop- 
erties, but we have included this also within the context 
of an impurity band model for comparison with the other 
cases. 
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TABLE I. Sample sizes considered for different values of x 
and p. 



We carried out simulations on lattices of linear size be- 
tween L — 7 and L — 14, which contained between 40 
and 110 Mn spins and between 4 and 21 carriers. The 
sizes that were considered and their labels are tabulated 
in Table || We averaged up to 700 samples per data 
point, depending on the size and temperature. Typi- 
cally at least 30-40 samples were averaged for each data 
point. In this section we present our results for global 
quantities such as the magnetization, the Binder cumu- 
lant (to determine T c ) and the magnetic susceptibility. 
We also present our results comparing ordered and dis- 
ordered samples. We focus on local quantities in the next 
section. 



A. Magnetization 

In Fig. we plot the average magnetization per Mn 
spin, SMn(T), as a function of temperature for all four 
combinations of x and p using samples containing be- 
tween 60 and 70 Mn spins in each case (thus finite size 
effects are similar in all four curves). Clearly, the critical 
temperature, T c increases with increasing x and p. The 
magnetization curves for cases B2 and C2, which have 
the same hole concentration px have similar numerical 
values, although the curve for x — 0.01, p — 0.3 appears 
to have a lower T c . The most important feature is the 
unusual shape of the magnetization curves. The magneti- 
zation decreases rapidly from full polarization at T = 0, 
leading to linear or concave upward shapes, ^similar to 
those found in the mean—field approximationQ'tj and for 
insulating II- VI DMSjHj'EJ and very unlike the strongly 
convex upward magnetization curves seen in conventional 
ferromagnets such as ironO 




FIG. 3. Mn magnetization as a function of temperature 
(normalized to saturation) for case A2 (x — 0.01, p = 0.1), 
case B2 (a; = 0.01, p = 0.3), case C2 (a; = 0.03, p = 0.1), and 
case D2 (a; = 0.03, p = 0.3). 

The carrier magnetizations, s c (T) for the same sizes 
shown in Fig. |^ are shown in Fig. The curves mirror 
the feature of the Mn magnetizations that increasing x 
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and p lead to greater polarization at higher temperatures 
- there is also a clearer distinction between curves for B2 
and C2, where x = 0.01, p = 0.3 polarizes at noticeably 
lower temperatures than x = 0.03, p = 0.1. One major 
difference from the results of the mean field study is that 
the curves for the carrier magnetizations appear to be 
much more like the Mn-magnetization curves, whereas in 
the mean field studju'tl the carrier curves remained al- 
most fully polarized until T was close to T c . This feature 
is in closer agreement with experiments than the mean 
field results. Ell 




FIG. 4. Carrier magnetization as a function of temperature 
for case A2 (x = 0.01, p = 0.1), case B2 (a; = 0.01, p = 0.3), 
case C2 (x = 0.03, p = 0.1), and case D2 (x = 0.03, p = 0.3). 



B. Critical Temperature 

We use the Binder cumulant curves, G(L,T) (see 
Eq. ([It])) to identify the critical temperature T c for each 
carr ier and Mn concentration. As described in Sec. 
|ll C| , T c is indicated by the simultaneous intersection of 
G(L, T) curves for sufficiently large sizes, L. Fig. || shows 
G(L, T) curve for x = 0.03 and p = 0.3 for both sample 
sizes (Dl and D2 of Table ^) over a wide temperature 
range (T/J = 0.02 to 0.6). As expected G(L,T) de- 
creases with increasing size L at high T whereas the re- 
verse is true at low T. At the transition temperature T c , 
G{L,T) is expected to be independent of L, which is in- 
dicated by the crossing point of the two curves, implying 
T c /J — 0.45. The solid curves are spline fits to the data 
appropriately weighted by the error bars. 

Figs. H - H show the G(L, T) data for the other three 
concentrations studied, each in the vicinity of the cross- 
ing point. Several hundred samples were generally aver- 
aged for each data point as indicated in the figure cap- 
tions. The solid curves are spline jits to the data over 
the entire temperature range studied (typically T/T c = 
to 3). For the one case where we studied more than two 
system sizes (x = 0.01, p = 0.1), which we exhibit in 
Fig. ^, all curves are consistent with a single intersection 
point, although a small size dependence cannot be ruled 



out .Ha Such a dependence would represent corerections to 
finite size scaling, arising from the relatively small sizes 
of the samples studied. (The effective linear size of the 

spin system is JV d 3 , which varies from 3.4 to 5). The small 
sizes, necessitated by the need to repeatedly diagonalizc 
the fermion Hamiltonian, also limited the dynamic range 
available in our study. Because of this, the curves for the 
sizes studied do not splay out very dramatically around 
T c . This, in turn, necessitates calculation of G(L, T) to 
high precision, which requires long runs and averaging 
over many samples. Despite these drawbacks, we find 
G(L,T) to be a more reliable estimator of T c than e.g. 
peaks in y(T) for the largest size studied. Based on Fig- 
ures H to H we estimate T c as shown in Table ||. 
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FIG. 5. Binder Cumulant as a function of temperature for 
x = 0.03, p = 0.3. The data shown is for cases Dl (Nd = 41, 
N h = 12) and D2 (N d = 61, N h = 18). 
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FIG. 6. Binder Cumulant as a function of temperature for 
x = 0.01, p = 0.1. The data shown is for cases Al (Nd = 53, 
N h = 5), A2 (N d = 69, N h = 7) and A3 (N d = 110, N h = 14). 
The curves for each size are obtained from a non-linear curve 
fit using the entire temperature range. 
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FIG. 7. Binder Cumulant as a function of temperature for 
x = 0.03, p = 0.1. The data shown is for CI (N d = 41, 
N h = 4) and C2 (iV<j = 61, iV h = 6). The curves for each 
size are obtained from a non-linear curve fit using the entire 
temperature range. Up to 700 samples were averaged for the 
data shown. 
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FIG. 8. Binder Cumulant as a function of temperature for 
x = 0.01, p = 0.3. The data shown is for N d = 53, N h = 16 
and Nd — 69, Nh = 21. The curves for each size are obtained 
from a non-linear curve fit using the entire temperature range. 

In Tabic [il] we compare the Curie temperature T c as 
determined in the Monte Carlo simulations here with the 
results of the mean-field approximation using Langevin 
functions for classical vector spins. 



X 


P 


T C MC /J 


rpMC jrpMF 


0.01 


0.1 


0.037 ± 0.004 


0.14 


0.01 


0.3 


0.12 ± 0.04 


0.48 


0.03 


0.1 


0.16 ± 0.02 


0.37 


0.03 


0.3 


0.45 ± 0.03 


0.82 



TABLE II. Comparison of T c determined using 
Carlo (T C MC ) and mean field (T C MF ). 



Monte 



As expected, the mean-field approximation overesti- 
mates T c . While the reduction due to fluctuations is 
only about 20% for the largest x and p studied, it can be 
much more significant (a factor of 5 or more) at densities 
at and below the metal insulator transition density. (We 
mention in passing that results in the mean-field approx- 
imation for quantum spins using Brillouin functions are 
higher than for classical spins by 60 — 80%). 

Since there are many models for DMS which are solved 
at the mean field level, having an understanding of the 
behaviour of these models when fluctuations are consid- 
ered is very important if one wants to make quantitative 
fits to data. For a model where Mn ions interact with 
carriers in an unperturbed host band, a similar -reduc- 
tion in T c was found in Monte Carlo simulations. E3 One 
important difference between the results here and those 
found in Ref. |l3| however, is that here mean field theory 
becomes more accurate with increasing carrier concen- 
tration at a given x, whereas the opposite appears to be 
the case in that model. 

The Monte Carlo simulations show clearly a strong de- 
crease of T c at low carrier density (px), in agreement 
with experiment, and rectify the unphysically large T c 
obtained in the position dependent mean-field treatment 
in this limit. From Table || it appears clear that the 
main dependence of T c comes from the carrier density 
(px), with a relatively weaker dependence on (x/p). 



C. Magnetic Susceptibility 

The Mn and carrier susceptibilities [Eqs. ( |i~5| ) and (fl£j)1 
are shown in Fig. ^ and Fig. [l^ respectively, as a function 
of temperature, for the same sizes and dopings that were 
used for Mn and carrier magnetizations in Fig. |^ and 
Fig.|. 




FIG. 9. Mn susceptibility per Mn spin as a function of 
temperature for for case A2 (x — 0.01, p = 0.1), case B2 
(x = 0.01, p = 0.3), case C2 (x = 0.03, p = 0.1), and case D2 
(a; = 0.03, p = 0.3). 
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FIG. 10. Carrier susceptibility per carrier as a function of 
temperature for for case A2 (x = 0.01, p = 0.1), case B2 
(x = 0.01, p = 0.3), case C2 (x = 0.03, p = 0.1), and case D2 
(as = 0.03, p = 0.3). 



At high temperatures, well above T c , one expects 



1 ^ 3k B T _ 
X ~ n(gfi B ) 2 S 2 \ T 



(26) 



for classical spins (the factor of 5 2 is replaced by 5(5 + 1) 
for quantum spins). We have investigated the behaviour 
of the inverse susceptibility in the temperature range 
T = 2 T c to 6 T c , and found that the effective spin 
ScS — 1.6 — 1.7 5, which suggests spin clusters (polarons) 
exist in the system well above T c . Fitting data close to 
T c (T c to 2 T c ) Eq. © gives 5 cff ~ 25, implying a sig- 
nificant enhancement of the Curie constant over the high 
T value. Ah. even larger enhancement has been seen in 
(Ga,Mn)N.y 



D. Effects of Disorder 

One of the major results of the mean field study was 
that samples with maximal disorder in the position of 
Mn spins wece found to have a higher T c than those with 
less disorder E This was explained as being due to carri- 
ers being able to lower their total energy more in regions 
with higher Mn density, by polarizing Mn spins and then 
hopping between these sites. However, it is expected that 
the mean field factorization, which assumes that the car- 
rier spin is either directed parallel or antiparallel to the 
overall magnetization, tends to align "islands" that might 
not otherwise be aligned until lower temperatures. This 
would be more likely to lead to a larger decrease in T c for 
disordered samples when thermal fluctuations are consid- 
ered, since global magnetization around T c is more likely 
to be destroyed than in ordered samples. 

We have tested whether the finding of disorder enhanc- 
ing T c is robust, by investigating two different cases. The 
first is for x = 0.01, p — 0.1, where we compared the 
weakly disordered case and fully disordered case, and the 



second is for x = 0.03, p — 0.3, where we compared the 
fully ordered case and the fully disordered case. 




0.04 0.06 

T/J 

FIG. 11. Binder cumulant for x ~ 0.01, p ~ 0.1 for the 
mildly disordered case, shown for Nd = 64, Nh = 7 and 
N d = 125, Nh = 13. 

In both cases we observed that the more ordered sam- 
ple had a lower T c . However, the enhancement is much 
smaller than obtained within the mean field treatment, 
as the arguments given above would suggest. 

For x = 0.01, p = 0.1, we determined the value of 
T c /J ~ 0.037 for the fully disordered case. We define 
mild disorder to correspond to the situation where Mn 
spins are chosen to be displaced from a fully ordered Mn 
lattice to one of the 12 nearest neighbour sites on the 
fee sublattice. We considered Nd = 64, Nh — 7 and 
Nd = 125, Nh — 13, and then used the Binder cumulant 
to estimate T c . 



b-b Disordered, N, = 61, N h = If 
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FIG. 12. Magnetization for x ~ 0.03, p ~ 0.3 for the com- 
pletely disordered case, shown for D2 (Nd = 61, Nh = 18) 
and the fully ordered case Nd = 64, Nh = 19. 
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Note that the number of Mn spins in our simulations 
are slightly different from those used for other cases, since 
an ordered lattice of Mn needs to be commensurate with 
the underlying lattice. Whilst we do not have a particu- 
larly accurate determination of T c , it is clear in Fig. ^ 
that the G(L) curves cross at a temperature well below 
the T c of the disordered case (T c /J = 0.037). 

In Fig. |lj we compare the magnetization calculated for 
x = 0.03, p = 0.3 with the completely random case D2 
and the completely ordered case with Nd = 64, Nh = 19. 
The x and p value were chosen since they give the largest 
T c and it was hoped that the effects would be more no- 
ticeable. It is clear that even when thermal fluctuations 
are included, the more ordered sample has a lower value 
of T c . (The magnetization at high T is a finite size effect) . 



V. RESULTS: LOCAL QUANTITIES 

We now discuss our results for various local quantities 
such as the local magnetization and local charge densities 
for various Mn and carrier concentrations. 



A. Charge densities 

In Fig. [l^ we show a histogram of the local carrier 
charge densities pi [see Eq. (fL8|)] at all Mn sites i, for 
x = 0.01, p = 0.1 (case A3) and averaged over 89 sam- 
ples at the temperature T/J = 0.01, compared with the 
corresponding distribution for x — 0.0093 andp = 0.1 ob- 
tained using the mean field approximation with Brillouin 
functions close to T — 0B Clearly the two distributions 
agree very well overall. The two peaks of the histogram 
suggest two populations of Mn sites. The peak at pi ~ 0.7 
is from sites which have a high probability of trapping a 
carrier, whilst the broader peak at much lower values of 
Pi is due to sites which have very little probability of hav- 
ing a carrier on them. For these sites, pi comes mainly 
through tailing from nea rby sites that have a high charge 
carrier density [see Eq. (fig)]. This shows that the charge 
carriers reside primarily on sites where there is a higher 
than average local Mn concentration, due to the strong 
interactions with these Mn spins. Because of the inho- 
mogeneous charge distribution, spins of Mn atoms on 
sites devoid of charge carriers have very small antiferro- 
magnetic couplings to the charge carriers, and therefore 
remain essentially free down to low temperatures. This 
explains the unusual shape of the magnetization curves. 

If we compare the histograms for the charge densities 
with the data obtained in the mean field model, we no- 
tice two important features. Firstly, the histogram is 
much less temperature dependent for T < T c . This is 
probably because we are generally looking at lower tem- 
perature scales than in the mean field case. Secondly, 
these distributions have large widths typical of highly 
disordered systems. This large width was first found in 



the mean-field case, and is the motivation for a simpli- 
fied phenomenologipal model for DMS, based on a two- 
component picture.E3 In this simplified model, one can 
divide the Mn spins into strongly and weakly interacting 
components depending on the temperature. This simpli- 
fied model has been shown to be adequate to reproduce 
the results of the full distribution at the mean field level, 
and can also explain experimental results on a qualitative 
basis BB 




o.oi o.i 

Local charge density (p.) 

FIG. 13. Histogram of local charge densities calculated us- 
ing Monte Carlo for x = 0.01, p = 0.1, N d = 110, N h = 11 
(solid line), and the corresponding distribution calculated us- 
ing mean field theory and Brillouin functions (Ref. 8). The 
temperature is T/J = 0.01 (T c /J ~ 0.037). 

In Fig. [l4| we show a histogram of the total charge 
densities p t at all Mn sites i, for x = 0.01, p = 0.3, 
Nd = 69, Nh — 21 and averaged over 41 samples at 
temperature T/J = 0.1. The most noticeable change as 
the temperature is varied is that the height of the peak 
at large pi decreases, corresponding to carriers becoming 
less localized at higher temperatures. 



0.03 




0.01 



0.01 

Local charge density (p.) 

FIG. 14. Histogram of local charge densities for x — 0.01, 
p = 0.3, N d = 69, N h = 21. The temperature is T/J = 0.1 
(T c /J ~ 0.12). 
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In Fig. |l5| and [l6| we show the histogram of the to- 
tal electron charge densities p, at each Mn site i, for 
x = 0.03, p = 0.1, Nd — 41, Nh = 4 at a temperature 
T/J = 0.14, averaged over 315 samples, and for x = 0.03, 
p = 0.3, N d = 41, N h = 12 at a temperature T/J = 0.4, 
averaged over 76 samples. (In both cases the tempera- 
tures are around 90% of T c ). Unlike the case of lower 
x, these distributions have only one peak, however sim- 
ilarly to the lower x case, the distribution is virtually 
independent of T below T c . 



3 0.02 - 
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FIG. 15. Histogram of local charge densities for x = 0.03, 
p = 0.1, N d = 41, N h = 4. The temperature is T/J = 0.14 
(T c /J ~ 0.16). 

Besides the change in the shape of the distribution, the 
width of the distribution is also significantly smaller at 
the higher hole concentration. This is consistent with the 
observation that the eigenstates at the Fermi energy are 
found to be delocalized at a = 0.03, whilst they appear 
to be localized at x = 0.01.1 



_r 0-03 






Local charge density (p.) 
FIG. 16. Histogram of local charge densities for x = 0.03, 



B. Local magnetizations 

We now consider the local magnetizations. In Fig. \\J\ 
we show a histogram of the local magnetizations M\° 
[see Eq. @] at all Mn sites i, for x = 0.01, p = 0.1, 
Nd = 110, Nh = 11 at three temperatures, one well be- 
low T c , one around T c /2 and one around 3T c /2. The 
number of samples to generate the histogram was 10 at 
T/J = 0.0053, 188 at T/J = 0.02 and 58 at T/J = 0.06. 

At the lowest temperature, there is a peak at around 
M iocai = q 7) and then 

a very broad tail that stretches 
to local magnetizations that are antiparallel to the over- 
all magnetization. For intermediate temperatures (T = 
0.02 J ~ Tc/2), there is evidence of two populations of 
Mn spins illustrated by the two peaks in the histogram, 
whilst for high temperatures (above T c ), there is a peak 
centered very close to zero local magnetization with some 
weight for small local magnetizations. This is consistent 
with the two component picture described in the sec- 
tion on the local charge distribution - at temperatures 
well below T c most of the Mn spins are strongly coupled, 
leading to a peak and then a tail of more weakly coupled 
spins. At intermediate temperatures, there are two com- 
parable populations of weakly and strongly coupled Mn 
spins, whilst at high temperatures, the local magnetiza- 
tions are all small and there is no long-range magnetic 
order. 




p = 0.3, N d = 41, N h 
(T c /J~ 0.45). 



12. The temperature is T/J = 0.4 



FIG. 17. Histogram of local magnetizations for x = 0.01, 
p = 0.1, N d = 110, N h = 11 at temperatures T/J = 0.0053, 
0.02, and 0.06 (T c /J ~ 0.04). 

Fig. [l8] shows a similar histogram of the local magne- 
tizations AfJ ocal , for x = 0.01, p = 0.3, N d = 69, N h = 21 
at temperatures T/J = 0.01, 0.033, 0.06, and 0.08 aver- 
aged over 12, 20, 70 and 106 samples respectively. The 
feature of two populations of Mn spins also appears to 
be present in this case, for T/J between 0.03 and 0.06 
(T c /J = 0.12). 
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FIG. 18. Histogram of local magnetizations for x = 0.01, 
p = 0.3, N d = 69, N h = 21 at temperatures T/J = 0.01, 
0.033, 0.06, and 0.08 (T c /J ~ 0.12). 




FIG. 20. Histogram of local magnetizations for x = 0.03, 
p = 0.3, N d = 41, N h = 12 at temperatures T/J = 0.04, 0.08, 



0.19, and 0.40 (T c /J ~ 0.45). 



Fig. |l9| and ^ show the temperature evolution of 
the corresponding histograms for x = 0.03, p = 0.1, 
(Nd = 41, N h = 4) and x = 0.03, p = 0.3 (N d = 41, 
Nh — 12) for temperatures ranging from well below to 
just below T c , For the former, with T c j J = 0.18 the data 
are for T/J = 0.01, 0.033, 0.08, and 0.14, averaged over 
49, 48, 48 and 334 samples respectively, whilst for the 
latter (T c /J = 0.45), we show data for T/J = 0.04, 0.08, 
0.19 and 0.40 averaged over 16, 22, 26 and 75 samples 
respectively. 
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FIG. 19. Histogram of local magnetizations for x = 0.03, 
p = 0.1, N d = 41, N h = 4 at temperatures T/J = 0.01, 0.033, 
0.08, and 0.14 (T c /J ~ 0.18). 

At each temperature there is typically a peak with 
some breadth, but unlike the two cases with x = 0.01, 
there is no double-peak structure seen. This can be un- 
derstood as increasing x leading to smaller relative den- 
sity variations for the Mn spins, and hence smaller fluc- 
tuations in the range of local environments, resulting in 
narrower distributions for local charge and magnetiza- 
tion. 



C. Correlation of charge and magnetization 

Large local magnetization of the local Mn spins are 
correlated with large local charge densities of the carriers 
at the Mn site. To illustrate this point, in Fig. [2l] we 
show the correlation of large local magnetizations and 
charge densities in a single sample with 110 Mn spins 
and 11 carriers at a temperature T = 0.01 J ~ 0.25 T c . 
Open circles correspond to sites with low charge density 
Pi < 0.1 and low magnetization M\° caX < 0.4. Solid filled 
circles correspond to sites with p t > 0.1 and M] ocal > 0.4 
Sites with < 0.1 but Nl\ ocal > 0.4 or with Pl > 0.1 



and M. 



local 



< 0.4 are shown as half-filled circles. Less 



than 14% of the sites are half-filled indicating a strong 
correlation between pi and Mj . 
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FIG. 21. Sites with large local magnetization and charge 
density in a sample with 110 Mn spins, 11 carriers at a tem- 
perature T = 0.01 J ~ 0.25 T c . 



Fig. f22| shows the joint distribution function for local 
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magnetization and local charge density for x = 0.01, p — 
0.1. The data inhabits a narrow band with the peak 
corresponding to localized carriers also corresponding to 
Mn spins with large local magnetizations. 
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FIG. 22. Joint distribution function of local magnetization 
and charge density for x = 0.01, p = 0.1, Nd = 110, Nh = 11 
at temperature T = 0.01 J ~ T c /4. 

In Fig. ^3] we plot the joint distribution function for lo- 
cal magnetization and local charge density for x = 0.03, 
p = 0.1. The distribution is similar to the previous case 
in that there is also a strong correlation between charge 
density and local magnetization. However, the data oc- 
cupy a smaller region of the Af] ocal , pi plane (which is 
not surprising since both distributions are narrower than 
in the insulating case). 
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FIG. 23. Joint distribution function of local magnetization 
and charge density for x = 0.03, p = 0.1, Nd = 41, Nh — 4 at 
temperature T = 0.06 J ~ T c /3. 



VI. DISCUSSION AND CONCLUSIONS 

In this study we have performed a Monte Carlo study 
of an impurity band model for III-V DMS. In order to 
do this we are simulating a model of fermions coupled to 
classical degrees (spins) and this requires diagonalization 
of the fermion problem for every classical spin config- 
uration. This is very time-consuming computationally. 
To speed up the procedure, we introduce quantum me- 
chanical perturbation theory coupled with Monte Carlo, 
which we call the PMC scheme. The method and tests 
are described in detail in Sec. Ill . Our model restricts 



charge carriers to an impurity band formed from the 
isolated acceptor impurity states introduced by the Mn 
ions, interacting via antiferromagnetic exchange with Mn 
spins, which we treat classically. This model is based on 
the picture for the low doping limit of localized carriers. 
Experimental evidence suggests that such an impurity 
band exists ip_.thjg vicinity of the metal-insulator tran- 
sition (MIT)E3'E3 The advantage of this type of model 
in comparison to other models that have been suggested 
for III-V DMS which start from a valence band point 
of view for the carriers, is that it naturally incorporates 
the physics associated with the MIT, which may be im- 
portant even in the metallic region. We have considered 
parameters appropriate for (Ga,Mn)As, although qual- 
itative features of the results may well apply for other 
III-V DMS. 

The impurity band model has previously been stud- 
ied using a mean-field approximationJZrtj in which it was 
found that the magnetic properties are very unlike those 
of a conventional ferromagnet. The magnetization curves 
Smu(T) were found to be concave upwards for a signifi- 
cant portion of the temperatures below T c , unlike conven- 
tional convex upward curves, and it was predicted that 
the magnetization should be inhomogeneous. The role of 
disorder in the material was also examined by comparing 
the magnetization curves for ordered arrays of dopants 
and impurities placed with varying degrees of random- 
ness. It was found that randomly placed impurities led 
to a higher T c than that found for an ordered lattice of 
impurities. However, in the mean field solution, the elec- 
tron spin is chosen to be either parallel or antiparallel 
to the total magnetization. Further, mean field approxi- 
mations neglect temporal fluctuations, and are known to 
overestimate T c . These overestimations could be signifi- 
cantly different in the ordered and disordered cases. 

The results of the Monte Carlo simulations confirm 
several of the results of the mean field study, while dif- 
fering on some others. Firstly, the shapes of the magne- 
tization versus temperature curves, Smu (T) are found to 
be unusual, as in the mean field study, but more linear in 
T compared with the strongly concave upward curves in 
the mean field approximation. Further, the carrier and 
Mn spin magnetization follow each other more closely, 
unlike the mean field, where the carriers remain polar- 
ized to much higher temperatures than the Mn spins. 



13 



We believe this difference is a direct consequence of the 
neglect of temporal fluctuations in the mean field approx- 
imation. The Monte Carlo data appear to be much more 
in accord with experimental results, which suggest that 
both Mn and carrier spins depolarize with roughly the 
same temperature dependenceEj In fact, the magnetiza- 
tion curves observed in this study bear a striking resem- 
blence to those obtained with magnetic circular dichro- 
ism results in Ref. [3l]. Whilst magnetic susceptibilities 
were not considered in the mean field study, in a Monte 
Carlo study applicable to II- VI semiconductors, a peak 
was found in the ordered phase, well below T c , in addi- 
tion to the singularity at T c .c9El Our susceptibility data 
show a single peak; however, the peak is at a tempera- 
ture significantly below the T c obtained from the G(L 1 T) 
curves. This is expected, since our sample sizes are small 
and the two peaks are not separated as a result. The 
relative proximity of the low T peak to T c compared to 
the II- VI case, makes this more difficult to resolve. Nev- 
ertheless, the explanation of the peak being due to free 
or partially-free spins that are outside the percolating 
magnetic cluster appears to be viable in both cases. 

To study the inhomogeneities in the magnetic be- 
haviour identified in the mean field study, we have stud- 
ied distributions of local quantities - charge density and 
magnetization. We have considered their joint distribu- 
tion as well, to investigate the correlation between the 
two. The calculations of the local charge density at 
low temperatures for the lowest density (x = 0.01 and 
p = 0.1), which appears to be insulating, are in very close 
accord with those obtained previously at the mean field 
levelia The generic feature that appears to be present 
is that the local charge density has significant disper- 
sion. For x — 0.01 we find a two-peaked structure in 
the local charge density distribution, corresponding to 
some sites having quasi-localized carriers and others hav- 
ing very low charge density, whilst for x = 0.03 there is 
still a broad distribution of charge densities, but there 
appears to be only one peak, rather than two, indicating 
that the carriers are in general delocalized. The distribu- 
tions of local magnetizations appear to reflect the same 
physical picture - for x — 0.01 at temperatures approxi- 
mately T c /2, there are two peaks, which appears to corre- 
spond to two populations of Mn spins - one which is still 
strongly magnetized, and one which is very slightly mag- 
netized, whereas for x = 0.03 this double-peaked struc- 
ture is not observed. The x = 0.01 behaviour was pre- 
viously predicted from mean field calculations,atJ and 
formed the basis for a phenomenological model that has 
been shown to have the capacity to describe experimen- 
tal magnetization curves SB Finally, the joint distribution 
of local charge density and local magnetization indicate 
that there is a strong correlation between higher charge 
density and large local magnetization. 

One of the main aims of the Monte Carlo study was to 
discover how important fluctuations are in determining 
the critical temperature T c at various x and p. As shown 
in Table In] mean field theory appears to be more precise 



in this model for increasing hole concentration, which is 
in contrast with Monte Carlo simulations on a valcnca 
band model where the opposite trend was observed lI§cZl 
A second item of interest is the quantification of the de- 
pendence of T c on the Mn concentration x and the carrier 
density px. For individual Mn coupled to free carriers, 
several models exist involving a combination of Fermi 
energy Ep ~ (px)3, and exchange energy E ex ~ Jx. 
The generic dependence of T c is found to be.oLiJie.form 
T c ~ x a (px)P, where mean field estimatesEjO'Ell give 
a = 1, with (3 varying between ^ and 1, whilst an analy- 
sis involving collective spin wave excitations with RKKY 
interactions yields a = 2, f3 = — ^ for weak coupling 
(Eex <C Ef) and a = — ^, f3 = 1 for strong coupling 
(E ex 3> Ep). That such a dependence exists for an im- 
purity band is not clear; nevertheless in the range of Mn 
and carrier concentrations studied (x = 0.01 - 0.03, px — 
0.001 - 0.009), our T c can be fit by an expression of the 
form T c ~ x a (pxY with a = 0.5 ± 0.15, (3 = 0.85 ± 0.15. 
As can be seen in Fig. |24|, a double logarithmic plot of 
T c /J versus x(px)^ yields a straight line with a slope 
of ^ . Since the dynamic range in carrier concentration is 
larger, the data restrict the range of allowed f3 more than 
a. It is interesting that both the dependence on carrier 
density for fixed Mn concentration (x), T c ~ p@ and on 
Mn concentration for fixed degree of compensation (p), 
T c ~ x Q+ ' 9 , yield exponents (3 — |, a + (3 — |, which 
lie within the range predicted by various treatments for 
free carriers, for which (3 varies from — ^ to 1 and a + (3 
from | to 2. Determination of the separate dependences 
on carrier and Mn concentration in experiment would be 
useful in clarifying the applicability of various theoretical 
models. 




FIG. 24. T c as a function of carrier concentration and Mn 
concentration for the four cases considered here. The line has 
slope | 

Experimental magnetization measurements also sug- 
gest that there is some degree of variation of local envi- 
ronment for Mn spins, as suggested here. When SQUID 
magnetometer measurements and anomalous Hall ef- 
fect measurements are compared, the magnetization 
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curves that are observed have very different character- 
istic behaviour. The SQUID measurements tend to be 
non-conventional magnetization curves as observed here, 
since they sample all Mn spins. However, the Hall effect 
measurements are made in transport and give curves that 
look quite conventional. It has been suggestedEd that this 
can be understood by the fact that the carriers that are 
measured in transport preferentially interact with the Mn 
spins that are strongly coupled to carriers, so that the 
magnetic properties of only one population of Mn spins 
are sampled in transport. The correlation between local 
charge density and magnetization seen here supports this 
scenario. 

There are still many experimental questions that re- 
main unresolved about the nature of the magnetic be- 
haviour in (Ga,Mn)As. A more systematic understand- 
ing of the dependence of T c on x and p will help deter- 
mine which models are more appropriate for which re- 
gions of the phase diagram. In particular, more accurate 
determinations of p appear to be one of the most im- 
portant ingredients. Recent steps in this direction where 
the carrier concentratioH-has been controlled with elec- 
tron doping are a start but this may introduce other 
complications due to dopant centers being spatially dis- 
tinct from the spins, unlike the case for Mn, where the 
two are in the same place. Local probes of the material 
like nuclear magnetic resonance will also help to uncover 
to what extent the magnetic environment is inhomogc- 
neous, and how this depends on p and x. 

Despite the qualitative agreements with experimental 
data, there are a number of effects iha,t are left out in the 
treatment of DMS by our modellrEII Within the mean 
field approximation, a number of hopping integrals have 
been considered-within a tight-binding description of the 
impurity bandja and it appears that the most impor- 
tant feature in determining T c is the density of states 
at the Fermi energy.Ej Other effects that have been left 
out in this model are carrier-carrier interactions, valence 
band states, spin-orbit effects and direct Mn-Mn interac- 
tions which could lead to frustration.E9 Whilst the direct 
Mn-Mn exchange should not be important for the x val- 
ues considered here, since there are relatively few Mn 
spins that are close enough for their antifcrromagnetic 
exchange to be important, they may become significantly 
more important at higher x > 0.1, and affect T c . 

In conclusion, we have performed Monte Carlo simula- 
tions on an impurity band model for III-V DMS. We have 
confirmed many of the features that were seen at the level 
of mean field simulations - unusual magnetization curves 
and inhomogeneous magnetization and charge density 
below T c . The unusual magnetization curves are also in 
qualitative agreement with many experiments. We have 
made a comparison of T c determined using each method 
and find that for larger values of Mn concentration and 
carrier concentration, the mean field determination of T c 
become more accurate. We also find a power-law relation 
between T c and carrier concentration and Mn concentra- 
tion, that could be compared with experiments. 
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APPENDIX A: EXACT SOLUTION FOR THE 
TOY MODEL 



Consider the Hamiltonian of Eq. ( j25|) 
j * 1 



ffCj/3, 



(Al) 



where 



N 



i=l 



is the total spin of the system. 

We parametrize S = 5(sin 9 cos (f>, sin 9 sin </>, cos 9) and 
perform the canonical transformation 



a! lT = cos -Cif + sin -e 40 c a , 

9 9 
dii = - sin — e l ^c n + cos -c a . 

It is then straightforward to show that 

JS N 



n 



As a result, we can diagonalize the Hamiltonian in the 
/c-space 



where 



E ka {S) 



-2t cos ka + '^jJoSa, 



(A2) 



(A3) 



with a = ±1 and the cyclic boundary conditions imply 
that ka = 2ttti/N, with n = 0, 1, ...,N — 1. There is a 
single lowest eigenvalue corresponding to k — 0, and then 
degenerate eigenvalues corresponding to left and right 
moving modes. 

Thus the grand-canonical partition function is [see Eq. 

®] 



N r r 1 
Z N = U dSU H(l + e-K E *"W-ri] 
„■ — i Lv J f„ _ 



(A4) 
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where we use the simplified notation J dVli = 

Jq smOidOi Jq 7 ' d(f>i for the integral over solid angle. One 
can avoid the 2N multiple integrals over individual spin 
angles, and replace them by an expression of the general 
form 



d(j> / d9sm6 I dSS 2 F N (S,e,(/)) 
ii Jo Jo 



-/3(£ fc „(S)-M)>| 



(A5) 



ho 



where 9, <f> define the orientation of the total spin S, and 
we use the fact that the magnitude of the tot al spi n, 5 , 
varies between and N. Comparing Eqs. ( A.4) and (A5), 
the definition of the "weight" Fn(S) is 



N j- , / N 

Fn(S) = J] / <Mi Sls-J2 s 



(A6) 



from which it is straightforward to derive the recurrence 
relation 

F N (S) = J <Kl N F W _ X (S - S N )9(N - 1 - |S - S w |), 

where the Heaviside function 9 insures that the argument 
of Fjf-i cannot have a magnitude larger than N — l. It is 
apparent that, in fact, Fjy(S ) = F^jS). This can easily 
be seen from the definition (A6) as well, since one can 
choose to define the angles O, with respect to the system 
of coordinates in which S = ST2 = Sz, and the result 
cannot depend on the particular orientation 9, <p of the 
total spin S. Using the new variable y — cos6*jv, and 
performing the trivial integral over ^>jv, the recurrence 
relation can be rewritten as: 



From the recurrence relation Eq. (A9), one can deter- 
mine the general solution for the weight function /n{%) 
for any integer N. This is a piecewise function, given by 



/at (a) 



2 N - 1 {N-2)\x 



\N-2 



k=0 



on the subinterval N — 2(m+l) < x < N — 2m of the sup- 
port interval [0, N], with m an integer < m < N/2 and 
where C% = N\/(k\(N — fe)!) is the appropriate binomial 
coefficient. 

The partition function is thus 



Z N = (4n) N dSS 2 f N (S)T\ (l + e 

J ° ka ^ 



-f3{E ka (S)-n) 



)■ 



where the chemical potential fi is fixed from the condition 
for the average number of fermions: 

(n c )= I dSS 2 f N (S)J2n g J-^ 

qa 







Zn 



where n qa — [e''(- B « o, ( s ') - W+l] -1 are the occupation num- 
bers of the fermionic levels. Since an analytical expres- 
sion is available for the weight function /jv(S'), these in- 
tegrals can be evaluated numerically. Once the chemical 
potential [i is known, any other expectation value, such 
as the total spin magnitude 

1 r N / \ 

{S) = Z^l dSS 3 f N (S)H(l + e-^-^) 



ka 



or the total fermionic spin 



iV 



Fn(S) — 2n J dy Fn -\(\J \ + S 2 — 2Sy) can be computed for any given temperature T. 



x9(N-l- Vl + S 2 - 2Sy). 

Defining 

F N (S) = (Anf^MS), 
we obtain the recurrence formula 



min {x+1 ,N — 1) 



\x-l\ 



dzzf N -i(z). 



This is supplemented by the "initial" condition 



1 



- Tx' 



(A7) 



(A8) 



(A9) 



(A10) 



which can be obtained through direct integration from 
Eq. @. 
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